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Abstract 

The impact of recently calculated next-to-next-to-leading order QCD corrections and soft 
gluon resummations on the extraction of higher twist contributions to the deep-inelastic 
structure function F2 is studied using the BCDMS and SLAC data. It is demonstrated 
to which extent the need for higher twist terms is diminishing due to these higher order 
effects in the kinematical region, 0.35 < X < 0.85 and > 1.2 GeV^ investigated. 
In addition, theoretical uncertainties in the extraction of higher twist contributions are 
discussed, and comparisons to results obtained previously are made. 



Deep-inelastic lepton-hadron scattering (DIS) data are an integral part of all global 
QCD analyses of parton distribution functions ||T], 0]. To guarantee that only the leading- 
twist part of the parton densities is extracted usually only data with sufficiently high 
momentum transfer and invariant mass W'^ = Q'^{l/x — 1) are selected, e.g., > 
2 GeV^ and W"^ > 10 GeV^ ||T|. These universal distributions can then be applied to 
other hard processes by virtue of the factorization theorem. Contributions of higher twist 
(HT) are expected to become increasingly important shown, e.g., in the 

infrared renormalon approach P], but they are suppressed by additional powers of l/Q^ 
with respect to the logarithmic dependence of the leading twist contribution. Hence 
HT operators can be safely neglected in all conventional global QCD analyses [jl], due 
to the aforementioned cuts on and, in particular, on W"^. 

However, twist-four contributions are interesting in themselves and can provide valu- 
able insight into higher quark-quark and quark-gluon correlators inside the nucleon. Not 
very much is presently know about these correlators apart from some studies in the frame- 
work of lattice QCD 0] . Recent measurements of the transversely polarized DIS structure 
function g2 may indicate that at least some HT are smallQ for an averaged of about 
2^3 GeV^ even at fairly large values of x. A better understanding of HT contributions, 
in particular their importance in describing low DIS data, is important in many re- 
spects. A wealth of low DIS data presently discarded in global QCD analyses would 
open up. Secondly, future experiments, e.g., the 'CEBAF @ 12 GeV program, focus 
on measurements at high x at comparatively low where HT are expected to be rel- 
evant. Finally, QCD analyses of longitudinally polarized DIS data would benefit. Here 
the available data are too scarce to allow for sufficiently 'safe' cuts on and W'^ when 
extracting leading twist polarized parton densities, and the importance of HT in relating 
the measured spin asymmetry Ai ~ gi/Fi (where HT may partly cancel in the ratio) to 
the structure function gi is still an open issue, see, e.g. Ref. 0. 

At present, data on the unpolarized structure function F2 obtained quite some time 
ago by the CERN-BCDMS and various SLAG experiments still offer the best testing 

^Measurements of 52 provide a unique possibility to disentangle its twist-three part since the leading 
twist part can be entirely expressed in terms of the rather well known structure function gi . 
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ground for HT studies since they cover a wide kinematical range up to large x and down 
to low values of with sufficiently good statistical accuracy. Therefore several attempts 
have been already made to disentangle leading and higher twist contributions to F2 [8-14]. 
As anticipated, HT effects become increasingly important as x — ^ 1. Based on a partial 
next-to-next-to-leading order (NNLO) QCD analysis it was argued [T^ that the need 
for HT is diminishing greatly when going from the NLO to the NNLO of QCD0. A first 



complete NNLO analysis |14| could not confirm this observation. However, in that paper 
no distinction between different sources of HT terms was made, and it could be that most 
of the remaining HT in NNLO is of purely kinematical origin and hence calculable (see 
discussion below Fig. 1). 

In the remainder of the paper we further extend these studies making use of the NNLO 



coefficient functions JT^ and recent estimates for the full x dependence of the NNLO 



splitting functions [|rB| based on the integer Mellin n moments calculated in |TP|. Another 



important ingredient of our analyses will be soft gluon resummations (SG) which have not 
been considered so far in extractions of HT contributions to DIS. The quark coefficient 
functions in DIS which link the quark distributions to the structure functions F12 
exhibit a large x enhancement of the form [ln'^^(l — x) /(I — where I = 1, . . . , 2/c in the 
0{a^) approximation to which needs to be resummed to all orders I^O], Finally, 



we compare the outcome of our analyses with results already available in the literature. It 
should be stressed that we mainly focus on a qualitative comparison of HT contributions 
to F2 extracted in NLO, NNLO, and NNLO including SG and do not attempt to provide 
a full set of parton densities or to extract as{Mz) from BCDMS and SLAG F2 data as 
was done, e.g., in [|, |10|. Instead we study the various sources of theoretical uncertainties 



in the determination of HT contributions to DIS. 

The DIS structure function F2 can be expanded in as 

= Ff (x, Q^) + (x, Q')/Q^ + 0{l/Q') (1) 
where Fj*'' denotes the contribution of twist-t. The leading twist part F^"* 'factorizes' 



^It should be noted that similar studies have been made also for the charged current structure function 
xF^ In it was also observed that NNLO contributions largely remove the need for HT at 

large x. 
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into a convolution of a perturbatively calculable coefficient function C2 and some non- 
perturbative parton density combination / which can only be determined from experiment 
so far: 

F!,'\x,Q^)=x[C,{a,{^^l),%,^)®f{^^),^^l)]{x) . (2) 

fij fit 

Once the non-perturbative input / is fixed at some reference scale fiQ its fif dependence is 
fully predicted by the well-known DGLAP evolution equation which schematically reads 

.2 ..2\ _ \^/„. I ..2\ ^ j:i ..2 ..2\ 



where /i/ and fir denote the factorization and renormalization scales, respectively. For 
simplicity we have limited ourselves in Eqs. (0) and @ to the non-singlet sector which is 
all what is needed for our analyses as will be discussed below. 

The coefficient function C2 in @ is known up to NNLO [l^ while only the ffist 



integer moments of the evolution kernels V in three-loop order have been calculated so 
far However, based on these results and further constraints on V estimates for the 



full X dependence of the NNLO kernels have been derived recently |[T8|. The residual 



uncertainties on the kernels were shown to be extremely small [0 in the large x region, 
X > 0.3, we are interested in. The relevant coefficients of the QCD beta function up to /?2 



which govern the running of as(/^r) are given in p2[. The appropriate matching conditions 



at flavor thresholds can be found in |2^. We are using the MS scheme throughout this 
paper. 

There are two sources of power corrections in contributing to F2{x,Q'^) in (|1]) 

beyond twist-2. The ffist one is of purely kinematical origin and can be entirely expressed 
in terms of the leading twist contribution F2'^\x,Q'^) in (^. It only takes into account 
effects of the so far neglected non-zero target mass M and approximately behaves like 
x^M^/Q^. Hence it gives a sizable contribution to F2 at large x whenever is of 0{M'^) 
or smaller. The 'target mass corrected' (TM) expression for F2^'*(x, Q^) is know for a long 



time and reads 24 
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with ^ = 2x/(l + ^/r) and r = 1 + Ax"^ M"^ / Q"^ . It is sometimes preferred to use F™ only 
after expanding in powers of M'^/Q'^ and retaining only the leading term. We will use 
both, the full and the expanded expressions for F™ in our fits. Any differences somehow 
represent part of the theoretical uncertainty in the extraction of higher twists from DIS 
data. 

The other source of power corrections cannot be related to F2^^(x, Q^) in general 
and provides important new insight into the QCD dynamics of higher quark-quark and 
quark-gluon correlations in the nucleon about which almost nothing is known yet. Also 
the evolution of these twist four operators has not been calculated yetQ. Therefore 
one has to fully rely on some model here. One interesting possibility is the infrared 
renormalon approach which allows to calculate explicitly the power corrections to the 
coefficient function C2 in i.e., to predict the x shape of F2^''(x, Q'^) up to some unknown 
normalization factor, which has to be determined from data, see p5| for details. This 
approach was used in some of the previous analyses of higher twists in DIS |TU| , |T^, [T^ . 



It turned out that, although the renormalon model could reproduce the general trend 
that HT increase as x — 1, the predicted x shape was somewhat off (see, e.g.. Fig. 2 in 
|T0|). Therefore we prefer not to fix the x shape of F2^'*(x, Q^) and instead use a simple 



multiplicative ansatz^ as was also employed in the first HT analysis by Milsztajn and 
Virchaux ^ as well as in |12, 1^. For each x bin of the data we simply introduce a factor 



Cht which is independent of Q^, 

Cut(x) 



F,{x,Q')=F™{x,Q') 



1 + 



Q 



2 



(5) 



i.e., F2^\x,Q'^) is approximated to have the same logarithmic dependence on as 
(a;,Q2). 

Finally, we also study the infiuence of soft gluon resummations for the coefficient func- 



Even if the evolution of twist four operators would become available, it would be extremely difficult 
to make use of it because the twist-4 operators which are accessible in DIS mix under evolution with 



other twist-4 operators which do not contribute to DIS and have to be taken from elsewhere 
;s F 



Sometimes F2'^\x,Q''') is defined in such a way that it does not exhibit any dependence, i.e., 
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tion C2 in (g) ||2y, ^ on the extraction of HT from DIS data. Since these resummations 
are operative also in the large x region where HT become important they may have a 
sizable impact on the size of the extracted HT coefficients Cht in ®. The SG take care 
of potentially large logarithms of the form [ln'~^(l — x) / (1 — where / = 1, . . . , 2/c in the 
0{a!l) approximation to C2 which threaten to spoil the convergence of the perturbative 
expansion by resumming them to all orders in a^. These resummations have been recently 



pushed up to next-to-next-to- leading logarithmic terms |21| and can be straightforwardly 
implemented in the Mellin n moment space which we also use for solving the evolution 
equations 

Before we turn to our numerical results let us specify the ansatz for the parton density 
function / in (|D and the data sets and cuts used in our analyses. Since we are only 
interested in the impact of NLO, NNLO, and SG corrections on the size of the higher 
twist coefficients Cht in (|^) we can limit ourselves to the large x region by taking into 
account only the BCDMS and SLAG data (using the BGDMS x binning) with x > 0.35 
101 . This is also the region where the estimates for the NNLO kernels V work extremely 



well |18|. As in [10-13] we then apply the non-singlet ('valence') approximation for F2 and 
only a single combination of parton densities / is required for proton target dataQ (we will 
comment below on a possible 'contamination' of the scaling violations of F2 at large x due 
to singlet and gluon contributions). In addition we select only data with > 1.2 GeV^ 
to stay away from the resonance region and the region where the perturbative expansion 
is expected to break down since becomes too large. Six bins in x remain after our cuts, 
so we have to determine six different HT parameters Cht (a; = 0.35), ... , CHT(a^ = 0.85) 
in (1^) from the fit to the data. 

We parametrize the input distribution f{x,fil) in (0) by the standard ansatz 

fix, fil) = Nx^il - x^il + 71 + 12X) (6) 

using fiQ = 1 GeV as the initial scale. It turns out, however, that all fits are rather 

insensitive to 71 and 72, so we will neglect these terms in the following. ^^(/Xr) is always 

^We refrain from using deuterium data as well. They are equally precise than proton data but may 
suffer from nuclear binding effects. We also discard SLAG data with x > 0.85 which are close to or 
in the resonance region. 
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computed by solving its renormalization group equation exactly in NLO or NNLO since 
the approximate formula is not sufficient for small scales. For simplicity we always 
choose /ir = /i/ = The evolution equations (§) are most conveniently solved in Mellin 
n space. There are, however, different ways how to actually truncate the solution at a 
given order k in perturbation theory (fc = 1: NLO, k = 2: NNLO), see, e.g., for 
a detailed discussion. On the one hand one can have a sort of iterative solution where 
all orders in still contribute, on the other hand one can impose strict power counting 
by keeping only terms up to Both approaches are widely in use in global QCD 

analyses and differ, of course, only by terms of 0{ai''^^^). This 'freedom' also represents 
part of the theoretical uncertainty and we will perform our fits using different solutions 
to (I). 

We should also mention that we have always added the statistical and systematic 
errors in quadrature as is commonly done in most QCD analyses of parton densities. In 
H] some of the systematic errors of the BCDMS data have been combined to a so called 



'main systematic error' |2^. The BCDMS data points are then allowed to float within 
this error with a common normalization factor to be determined from the fit The 
BCDMS data shifted in this way are used in some of the global QCD analyses. Unless 
stated otherwise, we refrain from using these data as they are strongly biased by the 
theoretical input that went into the analysis, and we stick to the original data sets |^ 
instead. In particular, it should be noticed that in ^ data with values as low as 
0.5 GeV^ were taken into account and that only the approximate solution for the running 
of as [0] was used which is known to be way off the exact solution at such small scales. In 



addition, it has been shown recently in ||T0[ that the introduction of the main systematic 
error is not a rigorous method to deal with systematic errors and that correlations should 
be fully taken into account to obtain reliable errors for the parameters extracted from a 
analysis. Nevertheless, as we are mainly interested in the qualitative effects of higher 
order corrections we treat all errors as uncorrelated which greatly simplifies the analysis 
and should not effect the general features of our results. The statistical meaning of the 
function is totally obscured in QCD analyses of parton densities anyway, and errors are 
very difficult to quantify ||30[. Therefore the errors for the Cht as obtained in our fits 
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Figure 1: Comparison of the results of our fits (solid lines) in NLO and NNLO of QCD 
and in NNLO including SG with the SLAC and BCDMS proton target data [0] (open and 
full circles, respectively) using «s(M^) = 0.117. The dashed lines refer to the target mass 
corrected part F™(x,Q'^) only with the additional higher twist factor (1 + Cnrix) /Q"^) in 

(2) 

Eq. (^ being omitted. The leading twist (LT) part F2 without TM is shown as well in the 
center part (dotted lines) for illustration. 



should not be taken too literally or be regarded as la-errors. As in [H, 12, 13| we allowed 



for a global normalization shift of the BCDMS data with respect to the SLAC data. 

In Fig. 1 we compare the results of our analyses of F2 according to Eq. (solid 
lines) with the SLAC and BCDMS proton target data [0] using a,(M^) = 0.117. Very 
good agreement with very similar values of is achieved for the fits in NLO and NNLO 
QCD as well as for a fit in NNLO including the SG mentioned above. The size of the 
required dynamical HT contributions can be inferred from comparing the dashed lines, 
which refer to the TM corrected part of F2 only with (1 + CHT(a^)/Q^) in Eq. (^) being 
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Figure 2: Extracted higher twist parameters Cut in six x-bins according to Eq. (H) for fixed 
and fitted values of as{Mz)- On the l.h.s. also the results obtained in [H are shown for 
comparison (crosses). For each x-bin the results are slightly horizontally displaced to avoid 
overlapping error bars. 



omitted, witli tlie full results (solid lines). Clearly the need for HT contributions is 
diminishing when going from NLO to NNLO. A further reduction, most noticeable at 
X = 0.85, is observed when the SG are taken into account as well. In the center part 
of Fig. 1 we also give the leading twist (LT) resuh Fg^'' without target mass correction. 
A comparison with the corresponding TM and full results shows that a major part of 
the HT contributions required to describe the data is of purely kinematical rather than 
dynamical origin and hence calculable. The extracted higher twist coefficients Cut for 
the six X bins are shown separately on the left hand side (l.h.s.) of Fig. 2. 

It should be mentioned that we have also tried to obtain as{Mz) from the fits rather 
than setting it to a fixed value. Unfortunately, the values we obtained turned out to be 
unacceptably small, as{Mz) ~ 0.105, way outside the current world average even within 
the error bars. The gross feature that the HT parameters get smaller when going to 
higher orders of perturbation theory is common to all results but less pronounced for 
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as{Mz) ~ 0.105 as can be seen by comparing the left and the right hand sides of Fig. 
2. This also demonstrates the strong correlation between the size of the Cht coefficients 
and the value of as(M^): the smaller the as{Mz) the larger the HT contribution as was 
already pointed out in ||T0[| . 



We have investigated the reason for not obtaining a more common value for as{Mz) 
from our fits in particular because reasonable values for as{Mz) have been obtained in 
previous analyses |T0|. The latter fits differ, however, in various aspects from our ones. 
In 1^ a NLO fit to all available BCDMS and SLAC proton and deuterium data was 
performed^ in the range > 0.5 GeV^ and 0.07 < x < 0.75. If we extend the range of 
data included in our fit to 0.07 < x, but still requiring > 1.2 GeV^, and introduce also 
a singlet and gluon contribution to F2, we obtain much more reasonable values for as{Mz) 
in NLO as well. This indicates that cannot be extracted from large x data alone and/or 
that singlet and gluon distributions are still important when studying scaling violations 
of F2 for X > 0.35, i.e., that a pure non-singlet approximation for F2 is not sufficient here. 
Another possible explanation for the small value of as{Mz) may come from our simplified 
treatment of the systematical errors. To investigate this further we have performed similar 
fits using the data in which include the BCDMS main systematic error shift mentioned 
above. In NLO this also leads to a significantly increased value for as{Mz) but almost no 
changes occur in NNLO. This does not come as a surprise since the main systematic error 
shift in 1^ was obtained only from a NLO fit, and it strongly depends on the theoretical 
input. It seems to be reasonable to assume that an as{Mz) of about 0.116 -^ 0.118 would 
be obtained in our fits once the correlations between different sources of systematical 



errors are properly taken into account as in ||10| , pj|] . This is however far beyond the scope 
of our more qualitative studies. Data at smaller values of x may also be important but a 
reliable NNLO analysis can only be performed once the complete NNLO kernels become 
available. As in |T2| , p!3[] , we therefore prefer to fix as{Mz) in all our fits to 0.117 rather 
than taking the fitted value which just seems to be an artifact of the approximations used 



^ As mentioned above, the approximate expression for as{nr) P7| was used in ||^ which is not vahd 
for smaU scales. In addition not all details of their analysis are specified in Therefore it is difficult to 
actually compare the result for as{Mz) (and the HT coefficients) obtained in [|| with our NLO results 
on a quantitative level. 
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in our analyses. 



Before we turn to a comparison with results obtained previously in the literature, let 
us briefly focus on the theoretical uncertainties in the extraction of Cht [and as{Mz)] 
from DIS data. We have already mentioned that there are different ways of how to treat 
the TM in (^) and that there is no unique way of actually solving the evolution equations 
(0) to a given order. We have performed several flts to take this ambiguity into account 
and it turns out that the gross features of CuTix) are basically unchanged but the precise 
values of the Cht (a;) do depend to a certain extent on the details of the analysis. For 
simplicity we have always chosen = f^f = Q but nothing prevents us from varying 
these two scale independently around Q which would also alter the values of the CHT(a;)- 
The dependence on the factorization scheme and on /ij can be removed if one expresses 
the scaling violations of F2 in terms of F2 itself using 'physical' evolution kernels, see, 
e.g., Ref. Another major uncertainty which is difficult to quantify stems from the 



form of the chosen ansatz (|^) for the HT terms (see the discussion above). In addition we 
cannot take into account the correct dependence of the HT terms since the relevant 
anomalous dimensions are not know (cf. footnote ^). Apart from studies of the scale 
dependence, see, e.g., Refs. |rT|, none of the other sources of theoretical uncertainties 



have been included in any of the combined extractions of HT terms and as{Mz) from DIS 
data. Since Cht (a;) and as{Mz) are highly correlated we therefore argue that errors on 
as{Mz) from those types of analyses are perhaps seriously underestimated and of limited 
use (in particular when one takes into account also the obscured statistical significance of 
in a QCD analysis |3^). To determine as{Mz) from DIS it seems to be much safer to 



introduce first appropriate cuts on and W'^ to remove the kinematical region which is 
contaminated by HT and then use the obtained as{Mz) value (or that from a global QCD 
analysis) to determine the CHT(a;) afterwards by relaxing the cuts and keeping as{Mz) 
fixed. 

Finally, let us compare our results with previous ones available in the literature. On 
the l.h.s. of Fig. 2 we also show the HT coefficients obtained in the NLO fit in p]. Despite 
the differences in both analyses (see also footnote P) reasonable agreement is achieved 
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with our NLO results for Cut- We cannot directly compare with the Cht in |T0|' [HI since 
their ansatz differs from (^), see also footnote ^, but the main features are the same, 
in particular the change of sign in Cht at around x ~ 0.45. As already mentioned in 
the introduction, it was argued in that NNLO corrections reduce the need for HT 
contributions to F2. However, in [|I^ only the NNLO coefficient functions were used and 
the parton densities were taken from a NLO global analysis without refitting them. Each 
X bin was assigned a separate normalization factor, called 'fioating factor', to substitute 



for the necessary refit of the parton densities [0], and the HT coefficients were taken from 



the IR renormalon model 11231 . Nevertheless, we confirm their conclusions. In addition. 



SG seem to further reduce the need for HT at large x. 

In Fig. 3 our fitted input distributions / are compared with the corresponding com- 
bination of parton densities from a recent global NLO QCD analysis which to a large 
extent excludes the kinematical region where HT become operative. Over the entire x 
range included in our analysis, 0.35 < x < 0.85, our NLO / differs from the MRST result 
by less than 10%. Only at NNLO with SG included, our / turns out to be significantly 
smaller than the MRST result as x — 1. This gives an indication of the importance of 
SG at large x. Since the resummations increase the coefficient functions smaller parton 
densities are needed in that region to describe the data. 

To summarize, we have studied the impact of higher order corrections on the extraction 
of HT contributions to the DIS structure function F2. The size of the required HT terms is 
diminishing when NNLO corrections are taken into account as was recently argued in the 
literature. However, even at NNLO with SG still a sizable, positive HT term is needed at 
large x while for x < 0.45 a small negative HT contribution is required. SG resummations 
tend to further reduce the size of the HT contributions and should be therefore included 
in future analysis of upcoming data at large x, e.g., from 'CEBAF @ 12 GeV. It was 
argued that simultaneous extractions of as{Mz) and HT terms in DIS may suffer from 
large theoretical uncertainties due to the strong model dependence and the unknown 
evolution of the HT contribution. 



11 




Figure 3: Comparison of our fitted input distribution / according to Eq. (|^) with the corre- 
sponding non-singlet (valence) combination of a recent MRST global analysis Shown are 
the deviations from the MRST 'central fit' at the input scale fif = fiQ = 1 GeV. 
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